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ABSTRACT 



A direct numerical simulation, based on spectral methods, has been used to 
investigate viscous, incompressible, steady, rotationally symmetric flow due to a sphere 
rotating with a constant angular velocity about a diameter. The equations of motion have 
been reduced to a set of three nonlinear second order partial differential equations in terms 
of the vorticity, the stream function and the azimuthal velocity. The calculations have been 
carried out for Reynolds numbers (Re) from the Stokes flow regime (low Re) to the 
boundary layer regime (high Re). 

The numerical results clearly show how the Stokes flow behavior for low Reynolds 
numbers, and the boundary layer behavior for high Reynolds numbers, are approached in 
the appropriate limits. Besides showing the flow streamlines, results have been presented 
for the torque and the skin friction behavior. It is shown that the present results are in 
excellent agreement with both available experimental data, and previously obtained 
numerical data. The radial equatorial jet which develops with increasing Reynolds numbers 
has been observed as expected from boundary layer collision behavior. No separation was 
observed for the range of Reynolds numbers considered, even near the equator. 
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I. INTRODUCTION 



The study of direct numerical simulation of steady viscous flow due to a sphere 
which rotates about a diameter with constant angular velocity in a fluid at rest has 
historically been the focus of a great deal of scientific attention. Most of the previous 
studies were based on finite-difference schemes. 

In this study, a different numerical scheme based on spectral methods was used. 
The validity of the numerical method and the results were checked with the previous 
studies. Spectral methods have become increasingly popular in recent years, especially 
since the development of fast transform methods, with applications in numerical weather 
prediction, numerical simulation of laminar and turbulent flows and other problems where 
high accuracy is desired for complicated results. 

The objective of this study is to compute the numerical simulation of a steady flow 
about a spinning sphere by using the spectral numerical methods. A Matlab program was 
developed by using spectral methods. However the method itself was not compared to the 
other methods used in previous studies. The reliability and validity of the numerical 
method were investigated to improve our current ability to compute and predict the 
evolution of flow for a particular geometry and conditions. 
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II. BACKGROUND STUDIES 



The spectral method is based on the collocation method which appears to have 
been used first by Slater and Kantorovic in specific applications in 1934. It was developed 
as a general method for solving ordinary differential equations. This method was revived 
by Wrigt in 1964. The applications of Chebyshev polynomial expansions to the initial 
value problems were involved in these studies. 

The spectral collocation method was applied to partial differential equations for 
spatially periodic problems by Kreiss and Oliger, who called it the Fourier method, and 
Orszag, who termed it “pseudospectral”. These were the earliest applications of spectral 
collocation or the pseudospectral method. This approach was very attractive because of 
its application to variable-coefficient and even non-linear problems. 

The Galerkin approach which depends on the same trial and test functions, was 
applied to a meteorological model by Silberman in 1954. This was the first serious 
application of spectral methods to PDFs. In 1970, Machenhauer and Rasmussen 
developed transform methods for evaluating convolution sums arising from quadratic non- 
linearity. Spectral Galerkin methods became practical for high resolution calculations of 
such non-linear problems. Applications in fluid dynamics were reviewed in the symposium 
proceedings edited by Voigt, Gottlieb and Hussaini in 1984. 

For the problem of the flow induced by a spinning sphere being considered in this 
study, the flow behavior for small values of the Reynolds number, Stokes [Ref. 1], Lamb 
[Ref 2], Bickley [Ref .3], Collins [Ref 4], Thomas and Walters [Ref 5], Ovseenko 
[Ref 6], and Takagi [Ref 7] have given theoretical results. For high Reynolds number, the 
boundary-layer equations have been studied by Howarth [Ref 8], Nigam [Ref 9], 
Stewartson [Ref 10], Fox [Ref 11], Banks [Refs. 12 & 13], Manohar [Ref 14], and 
Singh [Ref 15]. Over a large range of values of Reynolds number the flow has been 
investigated experimentally by Kobashi [Ref 16], Bowden and Lord [Ref 17], Kreith et 
al [Ref 18], and Sawatzki [Ref 19]. Another numerical method based on the use of 
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specialized techniques to obtain an approximation of the finite-difference form of the full 
partial differential equations was performed by Allen and Southwell [Ref. 20], Dennis 
[Refs. 21 & 22], Roscoe [Refs. 23 & 24], Spalding [Ref. 25] and Dennis, Ingham and 
Singh [Ref. 26]. 
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III. GOVERNING EQUATIONS 



A. DERIVATION OF EQUATIONS 



The governing equations for the problem are as follow; 



9u 

3t 



+ (u. V)u = -— Vp + vV‘u 



V.u = 0 



Define dimensionless variables; 






(3.1) Momentum equation 

(3.2) Continuity equation 




where a ; radius of the sphere, 

U : free stream velocity, 
p : density of the fluid, 

V : kinematic viscosity of the fluid. 

Introduce the dimensionless variables into Equation (3.1); 



Uv 9u* 



+ — (u*.V)u* 
a 



Divide Equation (3.3) by 



U.v 

a“ 



dn* 



U.a 



+ — 
V 



(u*.V)u* 



P a 



Vp* +v 



U 



V'u* 




+ vV^u* 



(3.3) 



(3.4) 



Hereafter the superscript * will be omitted for nondimension al terms. Define the 

U(2a) 

Reynolds number based on the diameter. Red = , and take the Curl of Equation 



(3.4) to be able to eliminate the pressure term and introduce vorticity, © = V x u where 
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Q)=lrCOt+leCOB+l0CGt:>. 



The transport equation for the <t>-component of the vorticity is 



d 2 f 

— co^i^-Vx(uxco^i^) = — I^V-- 



2 • 2 A 

r sin 0 



^ i (D 



(3.5) 



Convert the velocity in Equation (3.5) by defining the axisymmetric Stream function, V|/, 
and an “angular circulation”, Q. 



u = Vx 



¥ 

rsin6 



Q 

<1) ; Azimuthal direction (3.6) 

r sin D 



In spherical coordinates. Equation (3.6) is obtained as 





1 

QJ 

< 




1 8v|/' 


u = 


_r' sin 6 90 . 


’r- 


rsinS 8r _ 



'^■^rsine*'^ 



(3.7) 



where 



u = +Ueie+u^i^ 



(3.8) 



In the free stream vi/ = — Ur' sin' 6 or vi/ = — r' sin' 6 in nondimensional form. 

2 2 

By using the definition of the stream function it can be shown that 



(0:> = 



'a-vir 1 1 a 


^ 1 9\|/^ 


-1 'I 


dr' r" sin0 90 


k,sin0 90 ; 


Vrsin0 j 



(3.9) 
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Define operator D' as 



sind d f 1 d ^ 



to obtain the azimuthal component of the vorticity; 



Ob = - 



rsin9 



D^\|/ 



(3.10) 



(3.11) 



Call [ V X (u X Ciii>i<j>)].i(i, in Equation (3.5) nonlinear Term G(r,0,t) . By carrying out the 
calculations in spherical coordinates 



Vx(uxco^i^) 




a(\[/,D-\[t) 

3 (r,|i) 



+ 2D^\|/L(\|/) 



(3.12) 



is obtained, where 



|I = COS0 



L = 



ii a 1 a 

I-|i^ ar ^ r a|i 



a(z,,z2) az, az2 az, az2 

a(x,y) ax ay ay dx 



(3.13) 



Combine all terms to get 



a —2 Rcj 1 



a(\|/,D~\i/) 

a(r,)i) 



+ 2D^\|/L(\j/) 



= D> 



(3.14) 
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Equation (3.14) is called the Vorticity-Stream function form of the Momentum 
Equation in axisymmetric spherical coordinate system. For computational purposes, 
which will become clear later, define a modified stream function C(r,6,t), which is related 
to the usual Stokes stream function (\)/) by the relation; 

\|/ = rCsin0 (3.15) 

Also for any function f(r,9) , define operator D" such that 



where 



This gives 



where 



D‘(fr sin0) = rsin0D‘f 



D- = V- ^ 

^ 2 • 2 n * 

r sm 0 



co.=-D-C 



Q 



u = Vx(CiJ+^— i^ 
^ rsin0 



1 9 

u, =-^^3-(Csin0) 
rsm0 90 

1_^ 

r 9r 

Q 



Ue = — ~(rC) 



= 



rsin0 



(3.16) 

(3.17) 



(3.18) 



The modified stream function C is written as the sum C = C + c , where C = 0 , 
(which is the prescribed quiescent free stream condition), and c represents the disturbance 
produced by the presence of the spinning sphere. 
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Re 

Returning to Equation (3.5), rearrange by multiplying both sides with r^ ^ , 



, Rej a© , Re. , , 

-^G(r,0,t) + r'D^© 

r-D'C = -r'w 



(3.19) 

(3.20) 



where 



G(r,0,t) = [Vx(ux©^i^)]i^ 



(3.21) 



Note that the operator is defined as 



= 



_L_i 

r- 0r 



I' ar. 



+ ~ 



1 a r a ^ 



r" sin' 0 00 



1 



2 ' 2 f\ 

X Sin 9 



(3.22) 



Divide r'D^ in two parts which are and respectively and obtain 



r-D- =D? + 



— T7:De 

sin^ 0 ® 



(3.23) 



where Df and D^are as defined in Chapter IV. 

The nonlinear term in Equation (3.21), G(r,0,t), consists of 8 terms. These are 



1 a©,^ ac 


1 0©<j, 0C 


Cot0 8C 


Cot0 0©<„ 


r 0r 00 


r 00 0r 


“ r 0r 


■ r ^ 0r 


ac 


C 0©^ 


Cot0 0u^ 

2u ^ 


2u^ 0u^ 


r- 00 


r' 00 


r 0r 


r' 00 



Time integration of vorticity equation (3.19) is accomplished through the use of an 
explicit second-order Adams-Bashforth scheme for the nonlinear terms and an implicit 
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second-order Crank-Nicolson scheme for the viscous and linear terms. The calculations 
are made in physical r-space and spectral 0-space. If we denote vorticity (0(r,t) as the 
vector of N sine series then the discretized form of Equation (3.19) is 



.2 






- =f’^[«G,+PG,-.,]+[D;+A](^^ 



-I- CO 



I+Al 



At ■ 2 L-r 2 

where a and P are suitably chosen weighting parameters in the implicit scheme. 
Rearrange Equation (3.24) and obtain 



(3.24) 



- , Re. 




, , Re. 


D +A-r' / 


®,.A. =- 


D^+A-hr^— ^ 

. . 



CO 



.-r'Re,[aG,+pG,_„] (3.25) 



3 1 

where a= Y ,/? = -— and G = (v x (u x co<pi<j,)) - i^ 

[d^ -I-A]c(r,t + At) = -r'CO(r,t-i-At) 



(3.26) 



We find the governing equation for u<p in same way: 

3(4^, Q) 



dQ 1 

•-I- 



9t r" 8(r,Cos0) 






Re, 



where 



a 'T 

= Urt, and C = 



rSin0 



rSin0 



(3.27) 



and introduce D'Q = D'(u<„rSin0)= rSin0D‘u^ and divide each side by (r sin0) 
If Equation (3.27) is rearranged, the result is 



3u, 



8t 



3(ue>c) 

8(r,0) 



+ C^ 



— ^,Ln(rSin0) 



9(r,0) 



Re, 



D'u, 



(3.28) 
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The nonlinear term in Equation (3.28) is called H(r,9,t) and consists of 



1 9u^ ac 1 dC Cote ^ 

r3r30 re03r 



Cote 3u* dC 

~r^~dT " r' ae 



c 3u<i> 
r' ae 



uc) is written as the sum of two parts: potential and disturbance parts, 



_ ~ ■ - Sin6 

U<], — Uo +Ujj, where u<b = — 7- is the Stokes flow 

r 

solution which satisfies u* = 0 . The final form of Equation (3.28) is 



au^ 

at 



H(r,e,t)+-^D^(uJ 

KCd 



(3.29) 



ult,(t + At)-u]i,(t) 



2D^ 



At 



= ^[3H(t)-H(t-At)] 



Uo(0 + ui„(t + At) 



(3.30) 



ui,(t + At) 



'r^ Re, 


r'D-' 


= u„(t) 


r^Re, r-D-' 


r^Re, 


2At 


2 


2At ' 2 _ 


^ 4 



[3H(t)-H(t-At)] (3.31) 






At 



u^(t + At) = - 






At 



p2 Pp 

u.(l)-!-^[3H(t)-H(l-At)](3.32) 



1 1 



B. BOUNDARY CONDITIONS 



The boundary conditions for the numerical scheme are 



— 3c 


ac 




c = -C = 0 — 

3r 


= — r— = 0 u <i> = 0 at r= 1 

3r 


(3.33) 


c = 0 


CO = 0 U(|, = 0 at r = r„ 


(3.34) 



The Neumann Boundary Conditions in Equation (3.33) are handled with the 
Green’s Function Method that will be explained later. Zero boundary conditions in the 0 
direction are automatically satisfied with the choice of the sine expansions. 
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IV. METHOD OF SOLUTION 



A. NUMERICAL METHOD 

A pseudospectral or collocation method based on the technique developed by 
Marcus & Tuckerman [Ref. 29], is used to solve the equations of motion in time. Spectral 
methods may be viewed as an extreme development of weighted residual (MWR) which is 
generally used in a discretization scheme for differential methods. The trial functions and 
test functions are chosen as infinitely differentiable global functions. This is one of the 
features which distinguish spectral methods from finite-element and finite-difference 
methods. A pseudospectral transform method is used in this study. The approach taken in 
the transform method is to use the Fast Fourier Transform (FFT) to transform the 
functions to spectral space or Inverse Fast Fourier Transform (IFFT) to transform the 
equations to physical space. 

A physical process can be described either in time domain, by the values of the 
same h as a function of t, e.g., h(t), or else in frequency domain, where the process is 
specified by giving its amplitude, H, as a function of frequency f, that is H(f), with -<» < f < 
<». So once it is known whether the function is either odd or even, this information can be 
used in transformations. If h(t) is real and even then H(f) is real and even, if h(t) is real 
and odd then H(f) is imaginary and odd. It is necessary to go back and forth between 
these two domains by means of Fourier transform equations. As explained below, odd 
functions, using sine series expansions for 03 , c and u<j>, are used. The same properties are 
valid for derivatives of these expansions. Using these properties will increase 
computational efficiency. 

Functions are presented both in spectral space as a finite series of basis functions 
and by values at collocation grid points in physical space. The vorticity, 0), the modified 
stream function, C, and the velocity in the direction, u<^, are expanded as Chebyshev 
polynomials in the radial direction and as sine series in the transverse direction. The 
products and the derivatives in radial direction are done in physical space while the 
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derivatives in transverse direction are obtained in spectral space. Since each sine term in 
the expansion satisfies the homogenous 6 boundary conditions (to = 0, 'F = 0, u<j, =0 at 
6=0,7 t) exactly, no further 0 boundary conditions need to be applied. 

In general let a variable, say co, be expressed as 

N 

Q)(r,e,t) = 2cD„(z,t)Sin(6j (4.1a) 

n=l 

with 



M 

co„(z,t) = X(0„,„(t)T„(z) (4.1b) 

m=0 

where Tm(z) is Chebyshev Polynomial with -1 < z < 1 



0 = 

" N + 1 


n= 1,2,. 


..N 


(4.2) 


z = Cos 


'^mK^ 


m = 0, 1 ,. 


..M 


(4.3) 



An algebraic map is used to map the radial interval 1 < r < r„ to the interval 
-1 < z < 1 where 



r = l + L 



I + z 
b-z ’ 



lzl< 1 



(4.4) 



2L 

with b = 1 + (Note that to = 1 , which is the surface of the sphere.) 

r« -r„ 

L is a scaling parameter which is used to map the radial interval to the z-interval. 
A finite, but large outer r,., was chosen to avoid regularity problems in the radial 
differentiation. 

The sine expansions (4.1a) for co, C and u<t> satisfy the homogenous 0 boundary 
conditions and match exactly the periodicity and symmetry conditions in 0. The solution 
of the elliptic equation (3.20) is required to obtain the modified stream function, C, from 
vorticity coat any time level. Following [Ref 29], separable derivative operators, Dr^ and 
De‘, are defined as 
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where 



=Df + 7-D 

sin" 6 






De = 



(4.5) 



3r I 0r 






(4.6) 



In spectral 0-space the effect of operator De^ at any fixed radial location (physical r-space) 
can be written as 



D^f(6) = 



sin^ 0r-T + sin0.cos0 — - 1 



90 



90 



Ifj sin(j0) 



N 

...= If. 



J=> 



-(l-h-j^)sin(j0) + -j(l-hj)sin(20-hj0)+-j(l-j)sin(20-j0) 



(4.7) 



If this expression is written in 



the form of Dgf(0) = j^Syjfj 



f(i-2)(i-l) 



4 




(i + 2)(i + l) 



4 

0 



i - j-i-2 



i = j 



i = j-2 
otherwise 



(4.8) 



is obtained. 

This matrix is NxN if the representation for Dgf is truncated to N Fourier terms. 
Similarly, the operation of multiplying a function f(6) by sin'^(0) produces another 
NxN matrix which is called A matrix where the only non-zero terms are 



_ fal(i) = -i(i-i-l) 
■|a2(i) = -2i 



> = j 

• < j.i + j = even 



(4.9) 
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Time integration of the equations are obtained by using the scheme which is 
explained in Chapter III. 

The nonlinear terms, G and H, have to be carefully handled by transforming the 
spectral coefficients to physical space first, then multiplying the terms in physical space, 
and finally transforming the result of these products back to spectral space. 

The radial derivatives are evaluated by collocation methods and expressed as 
matrix operations on the vector of function values at the grid points in physical r-space. 
For (M+1) collocation points, we introduce the (M+l)x(M+l) matrix operation, that is 



D^'>(n) = D; + 





(4.10) 



where, I is the identity matrix. The operator [Dr^-r^(Re/At)+A] in (3.25) and (3.32) may 
now be written in block matrix form as 



D^'^(l) 

0 

0 

0 



0 

D^'^(2) 

0 

0 



a2(l)I 

0 

D^'\3) 

0 



0 

a2(2)I 

0 

D^”(4) 



(4.11) 



This matrix acts on the vector of a length of (M+l)xN for vorticity coefficients. 
Equations (3.25) and (3.26) are upper triangular, block matrix problems in physical r- 
space and spectral 0-space that can be solved by any solver with the Dirichlet and the 
Neumann boundary conditions discussed below. The equation (3.32) is treated in the 
same manner and leads to an upper triangular, block matrix problem. 

B. GREEN’S FUNCTION METHOD 

To enforce the radial boundary conditions in c-OJ equations, (3.25) and (3.26), 
Green’s Function Method has been developed which gives the correct boundary 
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conditions, (3.33) and (3.34), and allows the surface vorticity to develop naturally. If co 
and c consist of two parts, homogenous and particular solutions, the following is written: 



N 

® = and C 

j=i 



-Cp+X^jCj j=l,2,...,N 
j=i 



If (4.12) is introduced into (3.25) and (3.26), the homogenous parts as 
corresponding boundary conditions occur. 



E‘(Oj(r,0) = O (4.13) 



H^Cj =-r^o5j 



C0j(r = I,6i) = 5j. 
fOj(r = = 0 



€3(1 = 1,0.) = 0 

c.(r = oo,ei) = 0 



The particular parts with corresponding boundary conditions are 



E'( 0 „(r,e,t + At) = R(r,e,t) (4.15) H'c rr,6,t + At) = -r'(i3„(r,0,t) 



cOp(r = l,0)=O 



Cp(r=l,0) = -CU, 



ci)p(r = ~,0)=O 



Cp(r = oo,0) = O 



and are obtained from the block matrix defined in (4.11) 
To find A, 



is enforced. 



3c 3C 



(4.12) 
follow with 



(4.14) 



(4.16) 



(4.17) 
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If (4.12) is substituted into (4.17), 



^1 +T— I 

ar ^ dr ^=' ®' 






^1 

dr 



(4.18) 



is obtained. 

3c. 

= H. and -^U,e. = Aij (4.19) 

r=L0, 

N 

Finally from (4.18) and (4.19), H; = occurs and we solve for Xj as X = A*'H. 

j=i 

These X values can now be used in Equation (4.12) to arrive at the values of (oand c. 



Allow 



dC 
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V. RESULTS AND DISCUSSION 



Numerical calculations were performed for following parameters. 



Re 


M 


N 


L 


Too 


dt 


2 


64 


64 


7 


150 


0.05 


4 


64 


64 


6 


140 


0.05 


10 


64 


64 


4 


130 


0.03 


20 


64 


64 


4 


125 


0.05 


40 


64 


64 


4 


125 


0.05 


60 


64 


64 


4 


100 


0.03 


80 


64 


64 


3 


70 


0.03 


100 


64 


64 


2 


40 


0.01 


150 


64 


64 


2 


30 


0.03 


175 


64 


64 


2 


25 


0.03 


200 


64 


64 


2 


20 


0.05 


1000 


64 


64 


1 


17 


0.01 


2000 


64 


64 


0.6 


14 


0.01 


3000 


120 


128 


1.5 


15 


0.01 


5000 


128 


128 


1 


8 


0.01 



Table 1. Parameters Used in Numerical Calculations 

One of the main problems associated with many investigations that deal with 
complete numerical solution of Navier Stokes equations is the behavior of numerical 
scheme at high Reynolds numbers. It was found here that numerical results could be 
obtained for higher values of Reynolds numbers. For computational purposes, the results 
for Reynolds numbers up to 5000 are presented. 
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For comparison, T (Torque) is defined as 

T= Ja> ’SinGdS* 

S.r=l 



If (5.1) is nondimensionalized, the following occurs: 
= |o„Sin=ede , where o,. = 



_Ul> . 
9r r 



(5.1) 



(5.2) 



The small mapping parameter, L, with increasing Reynolds number provides a finer 
grid in the boundary layer. The re-circulation of the flow near the equator is another 
important result of a spinning sphere which will be explained later. Since most of the sign 
changes in velocity profiles occur in this region, decreasing L much below one is not a 
solution to gain some additional grid points in the boundary layer for high values of 
Reynolds number because finer grid near the sphere means coarser grid in the outer 
region. 

Several solutions for large Reynolds numbers have been obtained by integrating 
numerically the boundary layer equations found in [Ref. 13] which gives 



6.48 a'(0 

M = T 7 , where Re = (5.3) 

ReJ V 

There are some other investigations that recommend a different numerical factor 
such as 5.95 [Ref 8], 6.54 [Ref 12] and 6.53 [Ref 14]. In [Ref 26], Equation (5.3) was 
modified as 



M = 



6.48 



Re 



1/2 




(5.4) 



The Constants 6.48 and 3 1 in Equation (5.4) are for Reynolds numbers based on radius; 
these constants should be modified to 9.16 and 62 respectively for Reynolds numbers 
based on diameter as in this study. 
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Rcd 


Takagi 
[Ref. 7] 


Dennis, 
and Singh 
[Ref. 20] 


Banks 
[Ref. 13] 


Equation 

(5.4) 


Dennis, and 
Singh 
[Ref. 26] 


Present 


2 


50.307 


50.309 


6.48 


37.480 


50.305 


50.31 


4 


25.216 


25.218 


4.58 


20.082 


25.216 


25.23 


20 


5.401 


5.399 


2.05 


5.149 


5.398 


5.40 


40 




3.048 


1.45 


2.999 


3.048 


3.02 


100 




1.554 


0.916 


1.496 


1.554 


1.54 


200 




0.966 


0.648 


0.958 


0.966 


0.96 


1000 






0.290 


0.352 


0.348 


0.35 


2000 






0.205 


0.236 


0.234 


0.23 


3000 






0.167 


0.188 




0.19 


5000 






0.129 


0.142 




0.14 



Table 2. Nondimensional Torque (M) at Different Reynolds Numbers 
Figure 1 shows the variation of nondimensional torque with Reynolds numbers. It 
is seen that as the Reynolds number increases the general trend is quite consistent with the 
boundary layer solution given in [Ref. 13]. 




Figure 1. — , Banks [Ref. 13]; o, Sawatzki [Ref. 19]; +, Present 
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Figure 2 and Figure 3 show the transverse and azimuthal components of the 
dimensionless skin friction on the sphere, respectively. As it is in Figure 1, the higher 
Reynolds number, the closer it is to the boundary layer solution. Transverse and azimuthal 
components of dimensionless skin friction are defined as 



^^6 = 



^2 ' 




i^Re d , 


i UrJ 



(5.5) 



'^0 = 



2 TfduA 

V 






9r 



(5.6) 



^r=l 




Figure 2. Transverse Component of Skin Friction ; — , Banks [Ref. 13] 
o, Rea=20; +, Red=200; *, Red=3000 



22 




Figure 3. Azimuthal Component of Skin Friction; — , Banks [Ref. 13] 
o, Red=20; +, Red=200; *, Red=3000 

Present results agree well with the boundary layer solution near the pole but this is 
not so for the equator region. This is because the boundary layer solution is obtained by 
integrating parabolic partial differential equations from 0=0 with known initial conditions 
and the solution near the equator is determined from this procedure. Thus, it does not 
satisfy the physical boundary conditions. 

Since a finite r^ is used and the modified stream function is forced to be zero at r^ 
and on the sphere, a spurious re-circulating region of extremely weak flow (relative 
magnitude = lO’) is introduced to account for the continuity in the flow domain. The 
inflow at the pole changes to an outflow at the equator. The angular position of this 
change depends on Re but does not vary greatly in terms of radial distance from the 
sphere. However, the general trend is to get closer to the sphere with increasing Reynolds 
numbers. Figure 4 shows the stream lines in symmetric plane for different values of 
Reynolds numbers. 
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Figure 4. Stream Lines for Various Re 

In Figure 5, the jet effect can be easily seen at the equator for two different 
Reynolds numbers. The darker color represents the regions where the magnitude of 
dimensionless velocity in the symmetric plane is higher. The reference velocity in Figure 5 
is the velocity in the azimuthal direction on the sphere which has a value of unity. 




Figure 5. Jet Effect at Equator 
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Figure 6a and Figure 6b show the experimental flow visualization obtained by 
F. P. Bowden and R. G. Lord [Ref. 17], and present numerical solution for Red=3300 
respectively. Note that, unlike Figures, the velocity fields in Figure 6a and Figure 6b are in 
3-D. As it is seen there is a reasonable resemblance between experimental and numerically 
obtained pictures. In the experiment, the sphere was levitated by means of a magnetic 
field. It was noted in [Ref. 17] that, because of this magnetic field, there is a slight upward 
convection current so the smoke rises slowly; this explains why the top half of the picture 
is white and the bottom half black. Note also that, the length of the radial jet is relatively 
shorter in Figure 6a. This behavior can be explained by the presence of possible shear 
turbulence in the radial jet for high Reynolds numbers causing the laminar jet behavior to 
disintegrate. Since the present results are obtained by means of a numerical solution 
without any kind of turbulence modeling, there is no effect of turbulence in our cases, 
even for higher Reynolds numbers. 




Figure 6a. Smoke Picture of the Boundary Layer Flow for Re=3300. (From Ref. 17) 





In the boundary layer approximation, it is assumed that the velocity profile in the 
boundary layer just before the equator is the same as the one at the start of boundary 
collision. This assumption yields zero radial velocity at 0 = 7t/2, near the sphere. From 
the results it is seen that radial velocity is not zero at the equator near the sphere. Figure 7 
shows the variation of radial velocity with 9 for Re = 200 at different radial distances from 
the sphere. Note that the maximum value of radial velocity occurs exactly at the equator. 
Besides, there is an inward flow at both poles. Although the maximum value of radial 
velocity at the equator increases with Reynolds numbers, the thickness of the region where 
a general outward flow from the equator can be discussed decreases. 




Figure 7. Variation of Radial Velocity with 0 for Re = 200 
r= 1.0543; — , r=l.l 129; , r=1.1651 

Figure 8 shows the variation of radial velocity with radial distance at the equator 
for different values of Reynolds numbers. Details of this swirling Jet at the equator were 
given in [Ref. 27] and [Ref. 28]. 
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Figure 8. Variation of Radial Velocity with Radial Distance at 0=7t/2 o, Re=20 

+, Re=200; *, Re=3000 

Figure 9 shows the variation of dimensionless vorticity at the surface with 0 at 
different Reynolds numbers. The interest of this figure is to show that no separation 
occurs over the surface of the sphere, since a sign change must take place as a point of 
separation is passed. 




Figure 9. Variation of Dimensionless Vorticity with 0 
o, Re=20; +, Re=200; *, Re=3000 
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It is known that the boundary layer solution to be approached as Reynolds 
numbers increases, especially near the poles, and it is accurate to OCRe'"'^). Therefore, we 
- 1/2 

might expect that Re — at r = 1 and 6 = 0 will be of the form 



/ 



Re 



1/2 



9(0 ^ 



- 1/2 



(5.7) 



where A and B are constants. 



do) 



If Re' ■ at r = 1 and 6 = tc/ 2 is evaluted, it can be seen that this function is an 

OU 

increasing function of Re. It can be assumed that this expression is in the form of 

1 f 9(0^ 



Re ''H30 



,r=l 

^ 0=7C/2 



= CRe “ 



(5.8) 



where C and a are constants. 

The values of A, B, C and a were found by Dennis, Ingham and Singh [Ref. 26] to 
be 0.51, 0.56, 0.6 and Va respectively. The modified values of B and C for Reynolds 
numbers based on diameter are 0.79 and 0.5 respectively. 



R6d 


1/2 

Re,-"^(— ),., 

Ov 0=0 


0.79 

0.51+ 

Re/- 


1/2 

Re,-''^(— ),., 

C'V 0=jt/2 


0.5 Re,''" 


10 


0.695 


0.760 


0.43 


0.89 


20 


0.639 


0.687 


0.806 


1.06 


50 


0.597 


0.622 


1.15 


1.33 


100 


0.576 


0.589 


1.50 


1.58 


200 


0.552 


0.566 


1.79 


1.88 


1000 


0.528 


0.535 




2.81 


2000 


0.523 


0.528 




3.34 



Table 3. Calculated Values of Re''^’ 



dco 

le 



at r =1 , 0 = 0 and r = 1,6 = n/2 
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VI. CONCLUSION AND RECOMMENDATIONS 



In addition for the limiting cases of Re«l and Re»l for which Stokes flow and 
boundary layer solutions are available, the numerical results showed very good agreement 
with the those previously obtained by other numerical techniques. By changing some of 
the parameters such as r.. and L, the agreement with previous numerical results could be 
refined, although there is no proof that the other studies are correct within three decimal 
digits. 

Choosing the mapping parameter, L, should be done carefully so as to have many 
grid points in the regions of large gradients. A much more powerful mapping should be 
determined that provides a finer grid in two completely different regions, namely the 
boundary layer region and the re-circulation region (which though spurious, must be 
captured accurately by the numerical method for other calculations). 

A different FFT technique can be used that allows a number of grid points in the 9 
direction other than the ones that are a power of two. This would provide a relaxation in 
the choice of the number of grid points for different values of the Reynolds number which 
greatly affect the memory requirement. The problem is solved for 9 = 0 to 9 = K with the 
aid of sine expansions. Since the problem is also symmetric with respect to ti/ 2, by taking 
care of the even function properties it can only be solved for one quarter which simply 
doubles the grid points in 9 direction with the same memory requirement. 

The Green’s Function Method that deals with the Neumann boundary conditions 
should not be considered a trivial detail and must be used to have a stable convergence 
process. 

Obtained modified stream function values, C, can be easily used to solve the heat 
transfer problem due to convection from a spinning sphere. Derived additional energy 
equations are as follow 



29 



T{t + ^t) 



r^D^-Ec— 

At 



= -T{t) 



r-D~+Ec— 

At 



+ Ec^[ZK{t)-K{t-At)\ 



( 6 . 1 ) 



where Ec = Pr Re. 



The nonlinear term, K, in Equation (6.1) consists of 5 parts 

do dr do dr dO 



CrComf 

Since, for forced convection only. Equation (6.1) is uncoupled from Navier Stokes 
equations, the steady state values of modified stream function, C, can be used in the 
energy equation after solving the fluid problem which decreases the additional memory 
requirement due to different block matrixes in Equation (6.1). 
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APPENDIX A. PROGRAM STRUCTURE 



In the first part of the program, problem parameters were defined such as Re, dt, 
L, r_inf, M, N. Also calculations that are not dependent on time were performed in this 
part. To gain some additional memory Dl, D2 and D3 were calculated with D_Bmat 
function and saved as diagonal elements only. Non-diagonal elements of these matrixes 
were included in each iteration. 

Calculations showed that the matrix Ml (derivative operator in r direction) and 
consequently M3 and block matrixes that are diagonals of Dl, D2 and D3 are ill- 
conditioned. To avoid numerical errors in solving equation systems, the Singular Value 
Decomposition method was used. Necessary SVD decomposition of the matrixes and 
calculations of (O and c that are used in the Green’s Function Method were performed in 
this part. The built-in function in MATLAB, svd, was used for decomposition. 

Time dependent calculations basically consist of four steps. In the first step, Uo 
was calculated from Equation (3.32). In the second step, particular solutions of 0) and c 
were found from Equations (3.25) and (3.26). In the third step, H, which is used in the 
Green’s Function Method, was found and finally homogenous and the particular solutions 
were combined to get O) and c. 

Most of the calculations were performed in spectral domain which is a feature of 
the numerical technique. Only nonlinear operations were done in physical domain. 
Nonlinear terms were calculated with the functions nonlin and nonlinsp for co-c and u 
equations respectively. 

A convergence criterion was not used in the program. For high Reynolds Number 
(above 200), some over-shooting occurs before steady state. This situation might have 
caused wrong results if a criterion was used. Instead, the results were checked during 
calculations and the program was terminated after seeing steady state. 
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Solsvd : This function was used to solve equations after SVD decomposition. The 
limit used in this function reduces numerical errors due to ill-condition matrixes by 
skipping close to singular values in each iteration. 

Solsvd 1 : Since this function was used for solving homogenous parts in the 
Green’s Function Method and performed once before the time-loop, no limit was used in 
solsvd 1. 

Phy : This function was used for transformation of the variables from spectral 
domain to physical domain 

Operat : Derivatives in 9 direction in the spectral domain were performed with the 
function operat. Since these terms were used only in nonlinear terms, the output from this 
function is in the physical domain. 

Variables that were used in the program are 
L : Mapping parameter. 

Re : Reynolds Number. 

dt : Time interval between iterations. 

r_inf ; Outer boundary. 

step : Number of iterations. 

N ; Number of grid points in 0 direction. 

M : Number of grid points in r direction. 

xcor, ycor : Cartesian coordinates of grid points with respect to the center of the 
sphere (used for the presentation of the results). 

M 1 : Derivative operator in r direction. 

Dl, D2, D3 : Matrixes defined as in Equation (4.1 1). 
w : Vorticity. 

c ; Potential function as in Equation (3.15). 

u ; Velocity component in O direction. 

wtil, ctil : Homogenous part of the co and c defined as in the Green’s Function 
Method. 

Aij, H, landa : Variables defined as in the Green’s Function Method. 



32 



wp, cp : Particular solution of (oand c defined as in the Green’s Function Method, 
torq ; Dimensionless torque defined as in Equation (5.2). 
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APPENDIX B. PROGRAM CODES 



%% define program parameters 
clear all 

global EMI M2 M3 
L=7;Re=2;dt=0.05; 
r_inf= 1 50;step=4500; 

%% define grid points 



N=32;n=[0:N-l]'; 

teta=n.*pi./N; 

M=32;m=[0:M]'; 

z=cos(pi.*m./M); 

b=l+2*L/(r_inf-l); 

r= 1 +L. * ( 1 +z) ./(b-z) ; 

global N M 

xcor=cos(teta)*r'; 

ycor=sin(teta)*r'; 

E=emat(M); 
dispCE matrix is done') 

%% define Ml, M2,M3 

dz_dr=(b.*(L+r- l)-b.*(r- l)+L)./(L+r- 1 
B=diag(dz_dr); 

M1=B=^E; 

R=diag(r);global R Re dt 
M2=R*M1; 

M3=R*(2.*eye(M+ 1 )+M2)*M 1 ; 

%% define D1 D2 D3 
[D 1 ,D2,D3]=D_Bmat(M,N); 

% define initial w c wp cp u 
wp=zeros(M+ 1 ,N);cp=wp;w=wp;c=cp; 

dum l=sin(teta)*( 1 ./r.'^2)';dum2=[dum 1 ;zeros( 1 ,M+ 1 );- 1 .*flipud(dum 1 (2:N,:))]; 
dum3=fft(dum2); 
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u=wp; 



% define kmat (for teta der.) and cotteta (for nonlinear terms) 

dum=ones(l,M+l); 

kmat=n*dum; 

cotteta=[0;cot(teta(2 :N))] *dum; 
global kmat cotteta teta 

% svd decomposition for D1 D3 and define rsquare (for nonlinear terms) 

D 1 u=zeros(N*(M+ 1 ),M+ 1 ); 

D 1 1=D 1 u ;D 1 p=D 1 u;D3 1=D 1 u ;D3u=D 1 u ;D3p=D 1 u; 
rsquare=zeros(M+ 1 ,N); 
for i=l:N 

bas=(i-l)*(M+l)+l; 

son=bas+M; 

rsquare(:,i)=r.''2; 

Dl(bas,:)=[l zeros(l,M)]; % modify D1 for w @ r=inf 
Dl(son,:)=[zeros(l,M) 1]; % modify D1 forw @ r=l 
[D 1 l(bas:son,:),D 1 u(bas:son,:),D 1 p(bas;son,:)]=svd(D 1 (bas:son,:)); 
D3(bas,:)=[l zeros(l,M)]; % modify D3 fore @ r=inf 
D3(son,:)=[zeros(l,M) 1]; % modify D3 fore @ r=l 
[D31(bas:son,:),D3u(bas:son,;)<D3p(bas;son,:)]=svd(D3(bas:son,:)); 
end 

global rsquare 

% define c_ c_teta c_r u_ u_teta u_r 

c_=zeros(N,M+ 1 );c_teta=c_;c_r=c_; 
u_= 1 .*sin(teta)*( 1 ./r.''2)’; 
u_teta= 1 .*cos(teta)*( 1 ./r.''2)'; 
u_r=-2.*sin(teta)*(l ./r.'^3)'; 
global c_ c_teta c_r u_ u_teta u_r 

% define c boundary at r=l 

cbrl=zeros(l,N); 

% begin green function wtil ctil 

wtil=zeros(M+ 1 ,N''2);ctil=wtil; 
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for jj=l:N 
RR=zeros(M+ 1 ,N); 
duml=[zeros(l,jj-l) 1 zeros(l,N-jj)]; 
dum2=[duml 0 -l.*f]iplr(duml(2:N))]; 
dum3=fft(dum2); 
dum4=imag(dum3( 1 :N)); 

RR(M+ 1 , : )=du m4; 

wtil(:,jj*N)=solsvd 1 (D 1 1((N- 1 )*(M+ 1 )+ 1 :(N- 1 )*(M+ 1 )+ 1 +M,:), 
D 1 u((N- 1 )*(M+ 1 )+ 1 :(N- 1 )*(M+ 1 )+ 1 +M,:),... 
Dlp((N-l)*(M+l)+l;(N-l)*(M+l)^-l+M,:),... 

RR(:,N)); 



for i=N-l ;-l : 1 
sum=zeros(M+ 1,1); 
ii=(jj-l)*N+i; 
for j=i+l:N 

if rem(i+j,2)==2 

sum=sum+(-2*(i-l)).*wtil(:,(jj-l)*N+j); 

end 

sum(l)=0; 

sum(M+l)=0; 

end 

RR_=RR(:,i)-sum; 

wtil(:,ii)=solsvd 1 (D 1 l((i- 1 )*(M+1 )+ 1 :(i- 1 )*(M+ 1)+1+M,:),... 
Dlu((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Dlp((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 

RR_); 



tempwtil=(-l.*rsquare).*wtil(;,(jj-l)*N+l:jj*N); 
tempwtil(M+ 1 ,:)=zeros( 1 ,N); 
tempwtil( 1 ,:)=zeros( 1 ,N); 

ctil(;,jj *N)=solsvd 1 (D31((N- 1 )*(M+ 1 )+ 1 ;(N- 1 )*(M+ 1 )+ 1 +M,:)„ 
D3u((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,;),... 
D3p((N-l)*(M+l)+l;(N-l)*(M+l)+l+M,;),... 
tempwtil(:,N)); 



for i=N-l:-l:l 
ii=(jj-l)*N+i; 
sum=zeros(M+ 1,1); 
for j=i+l;N 

if rem(i+j,2)==2 
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sum=sum+(-2*(i- 1 )).*ctil(;,(jj- 1 )*N+j); 
end 

sum(l)=0; 

sum(M+l)=0; 

end 

wtil_=tempwtil(:,i)-sum; 

ctil(:,ii)=solsvd 1 (D31((i- 1 )*(M+1 )+ 1 ;(i- 1 )*(M+ 1 )+ 1 
D3u((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
D3p((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 

end 

end 

% find Aij matrix 



for jj=l:N 

blok=ctil(:,(jj-l)*N+l:jj*N); 
duml=phy(blok')'; 
dum2=(M 1 *dum 1 )'; 

dum3=[dum2; zeros(l,M+l) ; -l.*flipud(dum2(2;N,;))]; 
dum4=fft(dum3);dum5=imag(dum4( 1 ;N,:));dum5=dum5'; 
ctilr_ 1 (jj, 1 :N)=dum5(M+ 1 
end 

Aij=ctilr_l'; 

% find svd of Aij for landa solving 

Aij (!,:)=[! zeros(l.N-l)]; 

[Aiju,AiJs,Aijv]=svd(Aij); 

% begin time loop 

for ops=l:step 



if ops==l;prc=c;prw=w;pru=u;end 



% findRl 
Rl=zeros(M+l,N); 

R 1 (:,N)=D2((N- 1 )*(M+ 1 )+ 1 :(N- 1 )*(M+1 )+l +M,:)*w(:,N); 
for i=N-l:-l:l 
summ=zeros(M+ 1,1); 
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for j=i+l :N 

if rem(i+j,2)==2 
summ=sumin+(-2*(i- 1 )).*w(:,j); 
end 

end 

R 1 (:,i)=D2((i- 1 )*(M+ 1 )+ 1 :(i- 1 )*(M+ 1 )+ 1 +M,:)*w(;,l)+summ; 
end 

R1=-1.*R1; 

R 1 sp=zeros(M+ 1 ,N); 

R1 sp(:,N)=D2((N- 1 )*(M+ 1 )+ 1 :(N- 1 )*(M+1 )+ l+M,:)*u(:,N); 
for i=N-l;-l:l 
summ=zeros(M+ 1,1); 
for j=i+l:N 

if rem(i+j,2)==2 
summ=summ+(-2*(i- 1 )).*u(:,j); 
end 
end 

R 1 sp(:,i)=D2((i- 1 )*(M+ 1 )+ 1 :(i- 1 )*(M+ 1 )+ 1 +M,:)*u(:,i)+summ; 
end 

Rlsp=-l.*Rlsp; 

% find R2sp 

R2asp=nonlinsp(u,c); 

R2bsp=nonlinsp(pru,prc); 

R2sp=3/2.*R2asp-l/2.*R2bsp; 

R2sp=(- 1 *Re).*rsquare.*R2sp; 

Rtotsp=R 1 sp+R2sp; 

% find new u 

%% apply BC for u 

Rtotsp( 1 ,:)=zeros( 1 ,N); % BC for u=0 at r_inf 
Rtotsp(M+l ,:)=zeros( 1 ,N); % BC for u at r=l 

u(:,N)=solsvd(D 1 1((N- 1 )*(M+ 1 )+ 1 :(N- 1 )*(M+ 1 )+ 1 

Dlu((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 

D 1 p((N- 1 )*(M+ 1 )+ 1 ;(N- 1 )*(M+ 1 )+ 1 
Rtotsp(:,N)); 



for i=N-l:-l:l 
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sum=zeros(M+ 1,1); 
for j=i+l:N 

if rem(i+j,2)=2 
sum=sum+(-2*(i-l)).*u(:,j); 
end 

sum(l)=0; 

sum(M+l)=0; 

end 

Rtotsp_=Rtotsp(:,i)-sum; 

u(;,i)=solsvd(D 1 l((i- 1 )*(M+ 1 )+ 1 :(i- 1 )*(M+1)+1+M,:),... 
Dlu((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Dlp((i-l)*(M+l)+l;(i-l)*(M+l)+l+M,:),... 
Rtotsp_); 



end 



% findR2 

R2a=nonlin(w,c,u); 
R2b=nonlin(prw,prc,pru); 
R2=3/2.*R2a-l/2.*R2b; 
R2=(-l *Re).*rsquare.*R2; 
Rtot=R 1 +R2 ; 
prw=w;prc=c;pru=u; 



%% green function solution 
% find new wp 
%% apply BC for wp 

Rtot( 1 ,:)=zeros( 1 ,N); % BC for wp=0 at r_inf 
Rtot(M+l ,:)=zeros( 1 ,N); % BC for wp=0 at r_l 

wp(: ,N)=solsvd(D 1 1((N- 1 )*(M+ 1 )+ 1 ;(N- 1 )*(M+ 1 )+ 1 
Dlu((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
Dlp((N-l)*(M+l)+l:(N-l)*(M+l)+l+M,:),... 
Rtot(:,N)); 



for i=N-l:-l: 1 
sunn=zeros(M+ 1,1); 
for j=i+l:N 

if rem(i+j,2)=2 
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sum=sum+(-2*(i- 1 )).*wp(:,j); 
end 

sum(l)=0; 

sum(M+l)=0; 

end 

Rtot_=Rtot(:,i)-sum; 

wp(:,i)=solsvd(D 1 l((i- 1 )*(M+1 )+ 1 :(i- 1 )*(M+ 1 )+l+M,:),... 
Dlu((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:),... 
Dlp((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,;),... 
Rtot_); 



end 

% solve for new cp 
tempwp=-l ,*rsquare.*wp; 

tempwp(M+l,;)=cbrl ; % BC for cp=-c_r at r_l zero for pure spin 
tempwp(l,:)=zeros(l,N); % BC forcp=0 at r_inf 

cp(:,N)=solsvd(D31((N- 1 )*(M+ 1 )+ 1 :(N- 1 )*(M+ 1 )+ 1 
D3u((N- 1 )*(M+ 1 )+ 1 :(N- 1 )*(M+ 1 )+ 1 
D3p((N- 1 )*(M+ 1 )+ 1 :(N- 1 )*(M+ 1 )+ 1 
tempwp(:,N)); 



for i=N-l:-l:l 
sum=zeros(M+l,l); 
for j=i+l:N 

if rem(i+j,2)==2 
sum=sum+(-2*(i-l)).*cp(:,j); 
end 

sum(l)=0; 

sum(M+l)=0; 

end 

w_=tempwp(:,i)-sum; 

cp(;,i)=solsvd(D31((i- 1 )*(M+ 1 )+ 1 :(i- 1 )*(M+ 1 )+ 1 +M, 
D3u((i- 1 )*(M+ 1 )+ 1 :(i- 1 )=^(M+ 1 )+ 1 
D3p((i-l)*(M+l)+l;(i-l)*(M+l)+l+M,;),... 
w_); 

end 



% find landa 
% find H 
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duml=phy(cp’)'; 
dum2=(Ml*duml)'; % NxM+1 
dum3=- 1 .*c_r-dum2; 

dum4=[dum3;zeros(l,M+l);-l.*flipud(dum3(2:N,:))]; 

dum5=fft(dum4); 

dum6=imag(dum5( 1 :N,:)); 

dum6=dum6'; % M+lxN 

H=dum6(M+l,:);H=H'; % Nxl 
H(1)=0; 

landa=solsvd(Aiju,Aijs,Aijv,H); 

% find complete w & c 

dum I c=zeros(M+ 1 ,N);dum 1 w=dum 1 c; 
for jj= 1 :N 

blokc=ctil(:,(jj-l)*N+l:jj*N); 
blokw=wtil(:,(jj- 1 )*N+ 1 
dum2c=landa(jj).*blokc; 
dum2w=landa(jj).*blokw; 
dum 1 c=dum 1 c+dum2c; 
dum 1 w=dum 1 w+dum2w; 
end 

w=wp+dum 1 w; 
c=cp+dumlc; 

dum 1 =(phy(u’)+u_)'; 
dum2=Ml *dum 1 -dum 1 ; 
dum3=dum2(M+l,:)'.*(sin(teta).^2); 
torq ( ops )=-8 *p i/Re * trapz( teta, dum3 ) ; 
if ops>l 

err(ops)=torq(ops)-torq(ops- 1 ); 
end 

save tez4.mat u c w Re ops r_inf torq L N M dt; 

end % end of time loop 

function [Dl,D2,D3]=D_Bmat(M,N) 

global M3 R dt Re 

k=(M+l)*N; 

Dl=zeros(k,M+l); 
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D2=zeros(k,M+l); 

D3=zeros(k,M+l); 
for i=l:N 
for j=i;N 
if i=j 

Dl((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)... 

=M3+(- 1 *(i- 1 )*(i).*eye(M+ 1 )-Re/dt. *R.'^2); 

D2((i-l)*(M+l)+l:(i-l)*(M+l)+l+M,:)... 

=M3+(-l*(i-l)*(i).*eye(M+l)+Re/dt.*R.''2); 

D3((i-l)*(M+l)+I:(i-l)*(M+l)+]+M,:)... 

=M3+(- 1 *(i- 1 )*(i) *eye(M+ 1 )); 
end 
end 
end 

function E=emat(M) 

E=zeros(M,M); 

m=0:M; 

z=cos(pi.*m./M); 
for i=0:M 
for j=0:M 
if i~=j 

if ((i==0 I i==M) & (j==0 I j==M)) I ((i~=0 & i-=M) & (j-=0 & j-=M)) 
c=l;else 

if (i==0 1 i==M) & (j-=0 & j-=M) 
c=2;else 

if (i~=0 & i~=M) & (j==0 I j==M) 
c=l/2;end;end;end; 

E(i+ 1 J+ 1 )=c*(- 1 )^(i+j)/(z(i+ 1 )-z(j+ 1 )); 
else; 
if i==0 

E(i+lJ+l)=(2*M^2+l)/6; 
else 
if i==M 

E(i+lJ+l)=(2*M^2+l)/(-6); 

else 

E(i+ 1 J+ 1 )=- 1 *z(j + 1 )/(2*( 1 -z(j+ 1 r2)); 
end 
end 
end 
end 
end 
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function R2=nonlin(w,c,u) 

global kmat rsquare cotteta N M M 1 c_ c_teta c_r teta u_ u_r u_teta 
% variable list 

% 

% dw/dteta=non 1 dc/dteta=non2 dw/dr=non3 dc/dr=non4 (in phy) 
% du/dteta=non5 du/dr=non6 

% transpose matrixes for colomn operations in teta direction 

r=sqrt(rsquare)';w=w';c=c';u=u'; 

nonl=operat(w); 

non2=operat(c); 

non5=operat(u); 

wph=phy(w);cph=phy(c);uph=phy(u); 

non3=[Ml*wph']'; 

non4=[M 1 *cph']’;non4(: ,M+ 1 )=-0.5. *sin(teta); 

non6=[Ml*uph']'; 

term 1 =- 1 ./r. * non3 . *(non2+c_teta); 

term2=l ./r.*non 1 .*(non4+c_r); 

term3=-l.*cotteta./r.*wph.*(non4+c_r); 

term4=-l.*cotteta./r.*non3.*(c_+cph); 

term5=wph ./rsqu are' . * ( non2+c_teta) ; 

term6=(c_+cph)./rsquare'.*non 1 ; 

term7=l .*cotteta.*(non6+u_r).*(uph+u_).*2./r; 

term8=-l.*(non5+u_teta). *2. *(uph+u_) ./rsquare’; 

termO=term 1 +term2+term3+term4+term5+term6+term7+term8; 

ter=[termO;zeros( 1 ,M+ 1 ) ;- 1 . *flipud(termO(2:N,:))] ; 

RR2=fft(ter); 

R2=imag(RR2( 1 ;N,;));R2=R2'; 
function R2sp=nonlinsp(u,c) 

global kmat rsquare cotteta N M M 1 c_ c_teta c_r teta u_ u_r u_teta 
% variable list 

% 

% du/dteta=non 1 dc/dteta=non2 du/dr=non3 dc/dr=non4 (in phy) 

% transpose matrixes for column operations in teta direction 

r=sqrt( rsqu are)' ;u=u’;c=c’ ; 

nonl=operat(u); 

non2=operat(c); 

uph=phy(u);cph=phy(c); 

non3=[MPuph’]'; 

non4=[Ml*cph’]'; 

terml=-l./r.*(non3+u_r).*(non2+c_teta); 
term2= 1 ,/r.*(non l+u_teta).*(non4+c_r); 
term3=cotteta./r.*(uph+u_).*(non4+c_r); 
term4=-l.*cotteta./r.*(non3+u_r).*(c_+cph); 
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term5=- 1 .*(uph+u_)./rsquare'.*(non2+c_teta); 
term6=(c_+cph)./rsquare'.*(non 1 +u_teta); 
termO=term 1 +terrn2+term3+terrn4+tenTi5+term6; 
ter=[termO;zeros( 1 ,M+ 1 );- 1 .*fIipud(termO(2:N,:))]; 
RR2=fft(ter); 

R2=imag(RR2( 1 ;N,:));R2sp=R2'; 

function y=operat(x) 
global N M kmat 
re=zeros(N,M+l): 

X=re+i.*x; 

Y=i.*kmat.*X; 

Y_ful=[Y;zeros( 1 ,M+1 );flipud(Y(2:N,:))]; 

y=ifft(Y_ful); 

y=real(y(l:N,:)); 

function y=phy(x) 
global N M 
re=zeros(size(x)); 

X=re+i.*x; 

Y_ful=[X;zeros( 1 ,M+ 1 );conj(flipud(X(2:N,:)))]; 

dum=ifft(Y_ful); 

y=real(dum(l ;N,:)); 

function y=solsvd(u,s,v,b) 
w=diag(s); 

wmin=max(w)* 1 e- 1 2; 
for i=l:length(w) 

if w(i)<wmin;ww(i)=0;else;ww(i)=l/w(i);end 
end 

www=diag(ww); 

y=v*www*(u'*b); 

function y=solsvdl(u,s,v,b) 
w=diag(s); 
for i=l;length(w) 
ww(i)=l/w(i); 
end 

www=diag(ww); 

y=v*www*(u'*b); 
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